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COMPUTER  SIMULATION  OF  ELECTROMIGRATION 
IN  THIN  FILMS 

The  report  describes  an  on-going  effort  to  simulate  actions  which  physically 
occur  during  the  electromigration  process.  Important  characteristics  have  been 
included  such  as  diffusional  processes  which  account  for  lattice  mismatch  angle  or 
grain  boundary  energy,  grain  boundary  energy  equilibration  or  annealing,  and  the 
effects  of  an  encapsulating  layer.  This  description  is  meaningful  in  the  early  life 
of  films  when  inhomogeneities  of  temperature  and  current  density  are  not 


v 
.  v 


significant. 

The  author  points  out  the  formidable  nature  of  extending  this  simulation  to 
include  the  inhomogeneities  which  are  characteristic  of  films  near  end  of  life. 
Further  he  indicates  the  range  of  time  scales  involved  by  solving  for  "semi  steady 
state"  solutions,  because  of  the  shortness  of  diffusional  relaxation  times  as 
compared  to  film  lifetimes.  Future  work  in  the  area  of  Electromigration  physics 
include  improved  descriptions  of  film  growth  to  account  for  preferred  grain 
boundary  mismatch  angles,  and  investigations  of  the  quantum  mechanical 
interactions  of  conducting  electrons  with  point  and  extended  defects  such  as 


This  report  serves  a  double  purpose.  It  is  nominally  the  final  report 
for  the  renewal  project  under  USAF  contract  number  F-30602-81-C-0193  RADC  Task 
No.  N-6-5754  but  it  also  serves  to  review  the  progress  made  during  some 
fourteen  months  under  contract  SCEEE-PDP/85-0051  as  well  as  the  contract  just 
completed.  The  report  is  divided  into  four  sections.  The  first,  or 
introduction,  treats  quite  briefly  the  principal  background  based  on  earlier 
work  preceeding  the  SCEEE  contracts.  The  second  section  recounts  some  of  the 
complications  inherent  in  the  process  of  equilibrating  grain  boundary  tensions 
and  how  they  have  recently  been  resolved.  Section  three  deals  with 
electromigration  under  a  free  surface  and  section  four  is  concerned  with  the 
more  complex  situation  resulting  from  electromigration  under  a  passivating 
film.  Finally  we  discuss  briefly  what  the  succeeding  program  for  this 
investigation  might  be. 

I.  Review  of  Earlier  Work 

The  purpose  of  this  simulation  is  to  develop  a  representation  of  the 
actions  that  physically  occur  in  a  stripe  during  the  period  of  early 
deterioration  and  to  explore  the  influence  of  various  factors  in  stripe 
environment,  specifically  stripe  annealing  and  stripe  passivation  by  a 
constraining  surface  film. 

At  this  point  we  review  briefly  the  results  of  earlier  work,  preceding 
the  SCEEE  contracts.  The  first  step  was  to  develop  a  program  that  would 
generate  randomly  a  reasonable  network  of  grain  boundaries  for  a  stripe  of 
adjustible  dimensions.  For  this  an  extended  Voronoi  program  was  applied. 

Next  the  vertices  of  the  network  were  allowed  to  relax  to  attempt  to 
equilibrate  the  grain  boundary  stresses  as  it  was  believed  that  the  surface 
tensions  at  the  grain  boundaries  were  closely  linked  to  the  respective 
dif fusivities  and  that  their  equilibrium  at  the  vertices  would  tend  to  reduce 


Che  matter  build-up  or  reduction  with  time.  This  action,  as  a  form  of 
annealing,  should  tend  to  increase  stripe  lifetime. 

With  the  stripe  network  in  place  the  mass  motion  induced  by 
electromigration  was  introduced  through  the  appropriate  diffusion  equation 
modified  by  the  addition  of  a  driving  term.  It  was  apparent  from  the  values 
invoked  by  this  equation  that  the  time  constants  were  very  short  compared 
with  stripe  lifetimes  and  that  one  must  deal  almost  exclusively  with  the 
"semi-steady  stace"  solutions  in  calculating  the  mass  transfer.  The  cases  of 
the  free  surface  and  surface  constrained  by  passivating  layers  were  considered 
separately . 

1 I .  Finalizing  Stress  Equilibrium 

During  the  present  period  more  work  was  done  on  implementing  the  stress 
equilibrium  which  had  been  left  in  an  unfinished  state  earlier.  Although  the 
concept  of  allowing  the  network  vertices  to  move  to  relax  the  unbalanced  grain 
boundary  stresses  is  simple,  the  process  of  implementation  is  fraught  with  all 
the  difficulties  associated  with  finding  a  minimum  of  a  function  of  a  large 
number  of  coordinates.  In  this  case  it  was  not  necessary  to  set  all  net 
stresses  to  zero  but  only  reduce  them  to  below  a  certain  level.  Difficulty  in 
particular  occurred  when  stress  was  in  a  certain  grain  boundary  segment  strong 
enough  to  pull  together  the  segment's  end  vertices.  The  shortened  segment 
complicated  the  equilibrium  process  and  required  the  establishment  of 


carefully  considered  criteria  for  coalescence.  After  coalescence  there  was 


still  the  possibility  that  the  combined  vertex  would  separate  again  into  two 
vertices  creating  a  new  segment  roughly  orthogonal  to  the  original  segment; 
the  treatment  had  to  be  programmed  to  explore  this  possibility.  Finally  one 
had  to  allow  for  the  disappearance  of  any  triangular  grain  which  lost  a  side 
through  coalescence. 
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All  these  considerations  complicated  the  equilibrium  program  but  were 
eventually  mastered. 

III.  Mass  Motion ;  Free  Surface 

After  setting  up  the  grain  boundary  network  and  incorporating  the 
annealing  of  the  grain  boundary  tension,  the  action  of  forces  during  mass 
transport  was  investigated.  The  basic  equation  for  the  number  of  atoms  per 
unit  length  of  grain  boundaries  Mff(x,z,t),  in  stripe  numbered  a  is 


3M 
_ a 

3t 


A 


3M 

— 1 
a  3x  J 


+  D, 


(1) 


Here  D_l0  and  D| | 0  are  the  diffusion  constants  perpendicular  and  parallel  to 
the  dislocations  that  make  up  small  angle  grain  boundaries  in  the  network, 
i.e.  respectively  parallel  and  perpendicular  to  the  surface  of  the  substrate. 
The  quantity  Aa  relates  to  the  electromigration  drive  and  equals 
|  e  |  E  cos  <p  Z*/kT,  where  Z*  is  the  effective  charge  number,  a  dimensionless 
quantity  which  measures  the  strength  of  the  electromigration  force  and  <J>  is 
the  angle  which  the  grain  boundary  makes  with  the  current.  The  boundary 
conditions  for  the  solution  of  the  mass  flow  equation  is  that  the  flux,  given 
3M  (x) 

by  Dia(AaMa(x) - - - )  vanish  at  the  vertices  that  accumulate  material  at 

the  ends  of  the  segment  in  question.  On  the  other  hand  the  deficit  vertices 
(those  from  which  material  is  taken  away)  are  assumed  to  develop  voids,  which 
can  he  considered  most  simply  as  empty  circular  cylinders  going  from  surface 
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to  substrate.  The  boundary  condition  here  then  are  that  m(x,t)  be  zero  where 
we  take 


M(x,t)  =  M  +  m(x,t) 
o 


(2) 


Here  M0  is  the  equilibrium  number  "f  atoms  per  unit  length  in  the  grain 
boundary , 


M  =  N  d  6 
o 


(3) 


where  N  is  the  number  of  atoms  per  unit  volume,  d  is  the  thickness  of  the 
stripe  and  6  is  the  width  of  the  grain  boundary  taken  nominally  to  be  3x10“® 
cm.  Undoubtedly  the  concept  of  a  circularly  symmetric  void  is  an 
oversimplification.  At  one  time  we  considered  that  the  missing  material 
appeared  as  grooves  down  the  grain  boundaries  but  more  careful  consideration 
indicated  that  the  surface  driving  forces  would  strongly  favor  the  voids, 
although  they  are  probably  somewhat  irregular. 

Straight  numerical  solution  of  eq .  (1)  encounters  serious  convergence 
problems.  Recourse  to  analytical  approximations  soon  shows  that  the  natural 
time  constants  are  far  too  short  to  be  useful.  If  the  typical  segment  length 
is  L,  then  / D j_  is  of  the  order  of  an  hour  as  compared  to  stripe  lifetimes 

of  months.  Apparently  one  can  afford  to  disregard  the  transients  and  move  to 
what  we  have  cal  Led  "semi-steady  state"  solutions. 

In  this  section  the  case  of  deterioration  at  a  free  surface  will  be 
considered.  Excess  material  will  flow  unimpeded  up  through  the  grain 
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boundaries  and  spread  nearly  uniformly  over  the  surface  creating  hillocks  of 
low  contour.  The  matter  concentration  in  the  grain  boundaries  will  be 
strictly  time  independent  and  of  the  form, 

m(r)  =  m(x)  ( Tr/2d )  cos  (irz/2d)  (4) 

where  z  is  measured  from  the  substrate.  The  function  m(x)  satisfies  the 

equation , 


m 


(x)  -  Aq  m  (x)  =  (D | | /D±)  m(x), 


(5) 


Y1X  Y2X 

which  has  solutions  of  the  form  a^  e  +  a^  e  ,  where 


,  2  1/2 

Yl,2  =  1  1  Kd|  |/D1)(7r/2d)^  +  (§)  ] 


(6) 


However  it  turns  out  that  in  most  cases  A/2  <  (D^/D^)  (ir/2d)2  and,  in  the 
interest  of  simplifying  the  simulation,  the  A  has  been  dropped  or 


Y12=±  (Dli/Dl)1/2  (n/2d) 


(6a) 


Another  way  to  see  that  this  approximation  is  reasonable  to  look  at  the 
boundary  condition  at  the  vertices,  for  the  flux,  Dj^CAgMgCx)  -  m0'(x). 
Initially  the  M0(x)  is  a  constant  and  the  summation,  Fv  =  \  is 

positive  at  the  excess  vertices  and  negative  at  the  deficit  vertices.  Since 
m(v)«M0  it  is  reasonable  to  drop  the  term  Aa  m(x).  At  the  excess  vertices 
one  has  then  Fv  =  | IaDioma' (*)  I  • 

It  is  useful  to  designate  the  class  of  segment  by  the  nature  of  its 


bounding  vertices:  +  +  for  those  with  excess  vertices  at  each  end,  +  -  for 


those  with  an  excess  at  one  end  and  a  deficit  vertex  at  the  other,  and  -  -  for 
those  bounded  by  two  deficit  vertices.  For  the  +  +  segments  the  m(x)  is  given 
by  the  the  sum  of  two  exponentials.  For  the  +  -  segments  a  sinh  function 
appears  and  for  the  -  -  segments  m(x)  vanishes.  At  the  deficit  vertices  Fv 
gives  approximately  the  time  rate  of  increase  of  the  void  cylinders.  There 
is,  however,  a  small  correction  from  "spill-over"  from  the  +  -  segments. 

Since  the  sinh  has  a  non-vanishing  derivative  for  zero  argument,  there  will 
be  some  outward  diffusional  flow  (usually  small)  from  these  segments. 

There  is  a  universal  problem  in  all  this  simulation  in  determining  how  to 

divide  Fv  into  tne  respective  fa  that  flow  into  tne  individual  segments.  The 

exact  result  would  require  the  solution  of  a  large  set  of  interrelated  linear 

equations.  Rather  than  get  involved  in  this  complication  we  have  chosen  to 

3m  . 

move  by  successive  approximations.  Basically  f^  *  D  . j — - — |.  The  m^  are, 

11  OX 

however,  independent  of  i.  One  can  not  evaluate  the  derivative  without 
knowing  the  form  of  ma(x)  but  one  guesses  that  it  can  be  approximated  by 
mv/Ai.  Here  the  parameter  A^  is  set  equal  the  shorter  of  ,  the  length  of 
the  ith  segment,  and  the  parameter  (Dj^/D||)1/2  (2d/n).  Refinement  of  this 
approximation  can  proceed  by  iteration.  Next  one  determines  the  mv  as  equal 
to  Fv/XiDii/A,.  (  see  progress  report  #3  p.  4)  and  hence  m(x)  everywhere. 

The  end  result  of  these  simulations  can  be  presented  pictorially.  For 
this  we  show  the  network  at  a  particular  instant  in  time  with  scaled  circles 
around  the  deficit  vertices  and  graded  cross-hatching  on  the  +  +  and  +  - 
segments.  The  cross-hatching  represents  the  matter  concentration  in  the  grain 


boundaries  and  as  such  does  not  change  (after  a  brief  transient  period  which 
we  pass  over);  the  areas  of  the  circles  will  increase  linearly  with  time. 

This  situation  will  hold  throughout  the  period  of  initial  damage  and  will  only 
change  when  the  temperature  distribution  in  the  stripe  changes  because  of 
inhomogeneities  in  the  current  distribution  caused  by  the  presence  of  the 
voids . 


IV.  Mass  Motion;  Passivated  Surface 

In  this  section  we  consider  the  effect  of  the  presence  of  a  passivating 
layer  constraining  the  outflow  of  material  from  the  grain  boundaries. 
Mathematically  this  is  a  bit  more  complex,  since  now  the  concentration  of 
matter  In  the  grain  boundaries  is  building  up  with  time.  Again  we  disregard 
any  transients  and  look  for  "semi-steady  state”  solutions  of  the  mass  flow 
equation.  There  are  two:  m(x)  -  hx  and  m(x,t)  =  k  x2  +  At .  For  the  + 
-segments  only  the  first  is  applicable  corresponding  to  a  situation  where 
mac  ter  flows  uniformly  down  the  segment  with  no  change  in  concentration  with 
time.  For  the  ++  segments  both  forms  are  needed, 


m  (x)  =  k  x"  +  2D  ,  k  t  +  h  x. 
o  a  la  o  o 


We  take  the  Incoming  flux  at  the  left  and  right  vertices  of  the  segment  to  be 
given  by 


f.  =  D.  (k  L  +  h  , 
L/R  la  a  a  a 


where  the  origin  for  x  is  put  at  the  center  of  the  segment  of  length  L0. 
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!i$L 


Since  in 


=  k  (L  / 2)  +  2D  k  t  +  h  L  /2, 


it  follows  that 


ha  =  <mR  '  mL)/La 


k  =  (f,  +  f  J/2D,  L 
a  L  R  la  a 


A.  Initial  Approximate  Solution 


Again  one  faces  the  problem  of  avoiding  the  solution  of  a  large 


simultaneous  system  of  linear  equations.  We  have  chosen  to  take  eq.  (5)  and 


drop  the  h^  term  to  get  an  approximate  expression  for  k^  in  terms  of  m^, 


k  =  m  f(L/2)  +  2D,  t ] 

a  v  L  la  1 


By  substituting  the  approximate  expression  for  k^  one  into  eq.  (6)  obtains 


W»v  ■  L  -T-t^ - 

(T>  +2tDi0 


So  much  for  the  (  +  +)  segments.  Returning  now  to  the  (+  -)  segments  we 


have  (3rd  progress  report  p.  6) 


.  Vi  (II 

f  =  D,  m  /L  +  ~  L  ~ 
a  1  v  a  2  a  dt 


The  first  term  accounts  for  the  straight-through  flow  from  +  vertex  to  - 


vertex.  The  second  term  gives  the  mass  build-up  in  the  segment  that  comes 


with  the  steepening  of  the  slope  m^/L^.  Combining  we  have 


+  +  2 

F  =  y  f  =  y  D  .  L  f  (■=•)  +  2  D  t 

v  L  a  L  lo  aL  2  la 

a 


’  .  U  III 

+  \  {D  m  /L  +  3-  L  ) 

L-  lo  v  a  2  dt 

a 


The  evaluation  of  m  in  terms  of  the  known  F  is  straightforward  if  the  vertex 

v  v 

connects  only  (+  +)  segments  or  only  (+  -)  segments: 
a)  Only  (+  +)  segments 


+  +  L  -1-1 

m  =  F  [  l  D,  L  [(■=£)  +  2  D,  t]  ] 

v  v  L  L  io  o  1  2  lo  J  J 


(14a) 


which  indicates  an  m^  that  increases  steadily  with  time  as  long  as  all 
segments  stay  (+  +) . 


h)  Only  (+  -)  segments 


=  FV  L  I  <Dla/La>]  {1-exp-I  (2D,„t/Ln)/([V} 


io  a 


(14b) 


For  segments  of  modest  length  the  transient  will  probably  be  short  lived  and 

+  * 

m  soon  reaches  its  saturation  value  of  F f  ]  (D.  /L  )] 

V  L  °  X  O  0  J 

c)  The  equation  (13)  is  not  simple  to  integrate  if  both  (+  +)  and  (+  -) 

segments  are  involved  but  it  appears  that  again  m^  will  saturate  in  time  so 

that  can  be  considered  small  after  an  early  transient.  On  this  basis  we 
d  t 

take  m  =  m  -  m  „  where 
v  vl  v2 


m  =  F 
v  1  v  1 


K  D,„  <2D  1^  +  (T2)  )  +  l  (D.  /L  )]' 


,pr.  us  abbreviate  the  [  ]  by  n  (t).  Then 

m  =  F  h  (t) 
v  1  v  v 


(15a) 


-  ara^ 

and  F  =  m  tt  (t)  +  ( \  tr  L  )  -rr — 

\r  \t  u  /  n  n  r 


aim  a  —  in  «■  \  /  ■  \/  n  /  j 

v  v  v  c  2  a  at 

dm  ^mvj 

Since  the  final  term  is  small  we  approximate  =  — 


(13a) 


[iim 


t  d  7T 


k  Lc  Di»  t20!,' +  <r>  ] 


(16) 


o< 


„  +  -  L  +  +  L  L  2  .  -2 

Next  mv?  =  *—  [(  l  (/)  I  (/)(t  +  (/)  25“)  ] 
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Hopefully  (17)  is  a  small  correction  term. 


B.  Refinement  by  Iterations 


Next  one  pursues  an  iteration  procedure  to  refine  on  the  approximation  of 


dropping  h  at  eq.  11.  From  the  equations  (15)  and  (17)  one  knows  u»v  in  terms 


of  F  ,  L  ,  and  D,  .  From  the  mT  ,,  one  determines  h  from  eq.  (8)  and  new  (first) 
v  o  lo  L/R  a  n  v  / 


values  of  k  from  eq.  (10).  From  k  and  h  the  new  f'  can  be  obtained  for 
a  oo  vo 


the  (+  +)  segments  via  eq.  6.  For  the  (+  -)  segments  one  uses  eq.  12.  Since 

dm  „ 
v  Z 

-d<—  is  assumed  small,  it  is  probably  adequate  to  approximate  it  from  eq.  16. 


Next  one  applies  the  condition  [  f  =  F  ,  a  known  quantity.  Since 

o 


the  new  f^'s  do  not  automatically  satisfy  £fQ  =  Fy  and  are  indeed  propor¬ 


tional  to  m  ,  it  is  logical  to  renormalize  to  new  m  m'  . 

v  v  v 


HI  Til  '  r  + 

v  V  i  f  vo 


With  the  new  m"  one  is  ready  for  a  second  round  in  the  iteration  process. 


Brief  mention  is  here  made  of  the  changes  that  occur  at  the  deficit 


vertices.  The  <  0  continue  steadily  to  enlarge  the  voids  but  one  must  also 


look  to  the  spill-over  terms  which  tend  to  fill  them  up  as  discussed  in  the 


fourth  progress  report.  Here  the  spill-over  can  be  expressed 


—  + 

S  -  0  =  y  D,  trn/L 
L  lo  v+  o 
0 


Eventually  this  term  could  exceed  and  in  due  course  fill  up  the  void 
changing  the  vertex  from  -  to  +. 

Perhaps  the  simplest  way  to  check  on  the  progress  toward  convergence  in 
the  iteration  process  is  to  present  the  results  graphically  and  inspect  the 
superimposed  transparent  representations  of  successive  steps  in  the  iteration. 

C .  Time  steps 

After  determining  the  m  values  throughout  the  network  at  a  time 
sufficiently  short  to  avoid  non-linear  effects  from  saturation  and  spill-over, 
the  program  must  be  extended  through  subsequent  time  intervals  to  trace  the 
mass  motions.  Our  procedure  is  to  start  with  equation  (13a)  as  it  holds  at  two 
different  times  t^  and  =  t.  +  At.  Since  F  is  quite  time  independent  it 
drops  out  on  subtraction  and  we  have  our  relation  for  the  finite  difference 
of  m,  or  Am  =  mlt^)  -  m(t^): 

+  -  L  dm  dm 

-4mv  -  i«v(tl)S.  +  (  l  2^)((d~)t  -  <20) 

There  are  simplifications  for  two  special  cases: 

a)  If  the  +  vertex  links  only  to  -  vertices,  then  there  is  no  change  in 

m  unless  t,  lies  in  the  transient  region, 
v  1 

b)  If  the  +  vertex  links  only  to  +  vertices,  then  there  is  no  term 

d  m_ 

involving  - - . 

dt 


\  ^ 


3 


c)  If  both  +  and  -  vertices  are  involved,  we  proceed  by  considering  the 
derivative  terms  as  small, 


m^Ct^)  Ait 

*vl  wC  C  2 ) 


(21a) 


,  +  -  .  Am  m  (t  ) 

-Amv  »  -Amv  +  ir(t2>  \  l  ^-)  {-£•  -  (21b) 

Once  the  Am^  and  consequently  mv(t2)  are  determined,  one  can  proceed  to 

establish  the  new  h  and  k  and  to  calculate  the  additional  spill-over  in  At. 

a  a 

With  successive  time  steps  the  growth  of  m^  will  saturate  eventually, 
even  for  vertices  with  only  +  neighbors. 


V.  Continuing  Program 

Fortunately  we  have  some  support  from  an  industrial  company  to  continue 

somewhat  longer  our  investigation  of  this  simulation.  We  plan  to  implement 

with  a  program  calculation  the  analysis  for  the  passivated  surface  described 

in  part  4  and  to  compare  results  with  those  for  the  free  surface  using  both 

unannealed  and  equilibrated  networks.  Rupture  of  the  passivating  layer  at  m 

greater  than  a  certain  critical  value,  m  ,  should  also  be  considered. 

c 

We  plan  to  introduce  the  effects  of  surface  electromigration  in  inducing 


void  motion. 


In  the  long  term  an  effort  should  be  made  to  consider  the  influence  of 
temperature  and  current  inhomogeneities  caused  by  the  presence  of  the  voids 


although  this  is  a  formidable  complication. 
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